### This script will use a moving window analysis to try and differentiate between UL and Air data based on variance around the mean

####### First define the working directories

in.dir = "/home1/99/jc152199/varloop/"
setwd(in.dir)
out.dir = "/home1/99/jc152199/varloop/"

##### Import raw data, check the values, and examine the class of each column of the data frame

#air.data = read.csv('/home1/99/jc152199/MicroclimateStatisticalDownscale/ToAnalyse/MicroMacroMinMaxASCII.csv', header=T)
raw.data = read.csv('/home1/99/jc152199/underlog/underlograwdata.csv',sep=''), header=T)
head(raw.data)
str(raw.data)

##### Subset data to a few sites only for ease of processing

raw.sub = subset(raw.data, site=='WU9A2')
raw.sub2 = subset(raw.data, site=='WU11A2')
raw.data = rbind(raw.sub, raw.sub2)

##### Define function to perform calculations within loop and a list of sites to loop through

var.fun = function(x){return(sd(x,na.rm=T)*sd(x,na.rm=T))}
sites = unique(raw.data$site)

#Begin a loop that will calculate variance for a 24 hour moving window and plot that data

for (xsite in sites) {
sub.tdata = subset(raw.data, site==xsite)
sub.tdata$date = as.Date(paste(sub.tdata$year,sub.tdata$month,sub.tdata$day,sep="-"),"%Y-%m-%d")
start.date = min(sub.tdata$date)
end.date = max(sub.tdata$date)
sub.tdata$juldate = julian(sub.tdata$date, origin=start.date)
sub.tdata$datetime = strptime(paste(sub.tdata$year,"-",sub.tdata$month,'-',sub.tdata$day,' ',sub.tdata$time,':00',sep=''), format='%Y-%m-%d %H:%M:%S')
sub.tdata = sub.tdata[order(sub.tdata$datetime),]
vardata1 = NA
	
	if (nrow(sub.tdata)>1) {  
	pdf(file=paste(out.dir,xsite,".pdf", sep=""),width=11,height=7.5)
    par(mfrow=c(2,3), pin=c(3,3))      
    cat(xsite,'\n')    
		
		for (ii in 12:(nrow(sub.tdata)-11)) { 
		sub2.tdata = sub.tdata[(ii-11):(ii+11),]
		continuous = sd(sub2.tdata$juldate)
			
			if (continuous<0.55) {
			
			vardata = aggregate(sub2.tdata$logtemp, by=list(site=sub2.tdata$site), FUN = var.fun)
			vardata1 = rbind(vardata,vardata1)
			
			}
		
		names(vardata)[2]='variance'
		vardata = na.omit(vardata1)
		plot(sub3.tdata$day, sub3.tdata$maxdiff, cex.lab=.6, cex.main=.6, type='p', pch = 1, col=1, xlim = c(0,31), ylim = c(min(ranges),max(ranges)), xlab = "Day", ylab = "Microhabitat Temp Diff", main=paste(xsite," ",xyear,"-",xmonth,sep=""))
		points(sub3.tdata$day, sub3.tdata$mindiff, type='p', pch = 2, col=2)
		dev.off()
		
			
		}
	}                       
}
